{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysing Earth using EinsteinPy!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the eccentricity and speed at apehelion of Earth's orbit\n",
    "\n",
    "Various parameters of Earth's orbit considering Sun as schwarzschild body and solving geodesic equations are calculated\n",
    "\n",
    "### 1. Defining the initial parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "import numpy as np\n",
    "from einsteinpy.metric import Schwarzschild"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define position and velocity vectors in spherical coordinates\n",
    "# Earth is at perihelion\n",
    "M = 1.989e30 * u.kg  # mass of sun\n",
    "pos_vec = [147.09e6 * u.km, np.pi / 2 * u.rad, np.pi * u.rad]\n",
    "speed_at_perihelion = 30.29 * u.km / u.s\n",
    "omega = (u.rad * speed_at_perihelion) / pos_vec[0]\n",
    "vel_vec = [0 * u.km / u.s, 0 * u.rad / u.s, omega]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "   - Defining $\\lambda$ (or $\\tau$) for which to calculate trajectory\n",
    "     - $\\lambda$ is proper time and is approximately equal to coordinate time in non-relativistic limits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set lambda to complete an year.\n",
    "# Lambda is always specified in secs\n",
    "end_lambda = ((1 * u.year).to(u.s)).value\n",
    "# Choosing stepsize for ODE solver to be 5 minutes\n",
    "stepsize = ((5 * u.min).to(u.s)).value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Making a class instance to get the trajectory in cartesian coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/shreyas/Softwares/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/common.py:32: UserWarning: The following arguments have no effect for a chosen solver: `first_step`.\n",
      "  .format(\", \".join(\"`{}`\".format(x) for x in extraneous)))\n"
     ]
    }
   ],
   "source": [
    "starting_time = 0 * u.s\n",
    "obj = Schwarzschild.from_spherical(pos_vec, vel_vec, starting_time, M)\n",
    "ans = obj.calculate_trajectory(\n",
    "    end_lambda=end_lambda, OdeMethodKwargs={\"stepsize\": stepsize}, return_cartesian=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "   - Results are returned in SI units in cartesian coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Unit(\"s\"),\n",
       " Unit(\"m\"),\n",
       " Unit(\"m\"),\n",
       " Unit(\"m\"),\n",
       " Unit(dimensionless),\n",
       " Unit(\"m / s\"),\n",
       " Unit(\"m / s\"),\n",
       " Unit(\"m / s\")]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "obj.units_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "   - Return value is a tuple consisting of 2 numpy array\n",
    "     - First one contains list of $\\lambda$\n",
    "     - Seconds is array of shape (n,8) where each component is:\n",
    "       - t - coordinate time\n",
    "       - x - position in m\n",
    "       - y - position in m\n",
    "       - z - position in m\n",
    "       - dt/d$\\lambda$\n",
    "       - dx/d$\\lambda$\n",
    "       - dy/d$\\lambda$\n",
    "       - dz/d$\\lambda$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((13164,), (13164, 8))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ans[0].shape, ans[1].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating distance at apehelion\n",
    "   - Should be 152.10 million km"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$1.5205967 \\times 10^{8} \\; \\mathrm{km}$"
      ],
      "text/plain": [
       "<Quantity 1.52059668e+08 km>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = np.sqrt(np.square(ans[1][:, 1]) + np.square(ans[1][:, 2]))\n",
    "i = np.argmax(r)\n",
    "(r[i] * u.m).to(u.km)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Speed at perihelion should be 29.29 km/s and should be along y-axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$29.300051 \\; \\mathrm{\\frac{km}{s}}$"
      ],
      "text/plain": [
       "<Quantity 29.30005147 km / s>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((ans[1][i][6]) * u.m / u.s).to(u.km / u.s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating the eccentricity\n",
    "   - Should be 0.0167"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.01648174751185976"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xlist, ylist = ans[1][:, 1], ans[1][:, 2]\n",
    "i = np.argmax(ylist)\n",
    "x, y = xlist[i], ylist[i]\n",
    "eccentricity = x / (np.sqrt(x ** 2 + y ** 2))\n",
    "eccentricity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting the trajectory with time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# Changing Backend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/shreyas/Softwares/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/common.py:32: UserWarning: The following arguments have no effect for a chosen solver: `first_step`.\n",
      "  .format(\", \".join(\"`{}`\".format(x) for x in extraneous)))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fdecba49ef0>,\n",
       " <matplotlib.lines.Line2D at 0x7fdecba5d470>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAGBCAYAAACTjk3QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl0XOd55/nvUysK+75vBMEN3EmIErVLlm15iSTbsi3b7diO3Uq643R3ejJn3OMzSZ/MnDNJ58yZnu7ESRS3YyfteJdt2ZZl7bsoEdxJECQAYt/3HYVa3vmjijAoAQRAoupeVD2fc3BYy2Xd3y0A9eB97/u+V4wxKKWUUstxWB1AKaWUfWmRUEoptSItEkoppVakRUIppdSKtEgopZRakRYJpZRSK0qaIiEi3xKRQRE5v4Zt7xaRkyISFJFH3/XcMyIyLiK/jF1apZSyh6QpEsC3gQfXuG0n8EXgX5Z57q+Az29MJKWUsrekKRLGmFeB0aWPicjWaMvghIi8JiI7o9u2G2POAuFlXucFYCouoZVSymIuqwNY7AngD4wxzSJyK/AN4H6LMymllG0kbZEQkXTgduBHInL1Ya91iZRSyn6StkgQ6WobN8YcsDqIUkrZVdKck3g3Y8wk0CYinwSQiP0Wx1JKKVuRZFkFVkS+B9wL5AMDwJ8BLwJ/C5QAbuD7xpg/F5FbgJ8COcA80G+M2R19ndeAnUA6MAJ82Rjzm/gejVJKxUfSFAmllFLrl7TdTUoppVanRUIppdSKkmJ0U35+vqmurrY6hlJK2caJEyeGjTEFq22XFEWiurqahoYGq2MopZRtiEjHWrbT7iallFIr0iKhlFJqRVoklFJKrUiLhFJKqRVpkVBKKbUiLRJKKaVWpEVCKaXUirRIKKWUWpEWCaWUUiuyVZEQkW+JyKCInF/h+XtFZEJETke//jTeGZVSKpnYbVmObwN/DfzTdbZ5zRjz0fjEUUqp5GarImGMeVVEqq3OodTNmB2fIxQIEQ6GCYciXy6vi6ziTADGeiYwYYOI4HAJTrcTj8+NJ9VjcXKl3stWRWKNjorIGaAX+BNjzIXlNhKRx4HHASorK+MYTyWS4EKIuYk5QsEw2SWRD/mWN9oY75skMBcgMB8kMB8goyCdI48dBKDhR2eYGZ295nXyt+Ryy6cil1M/9bPz+Kf91zxfvKOAg4/sBeClb7wBgDfNgyfVgyfVTf6WXErrigGYHJzGl+nF5XUhIrE7eKXYfEXiJFBljJkWkQ8DPwO2LbehMeYJ4AmA+vp6vfyeWpExBv/0AikZXgDa3ulksGWYmbFZ/NMLAKRm+7jn948CMDU8w/yUH0+qm5TMFNwpLtLz0hZfb+f9tYQCYRwuBw6n4HA68C5pJRx8eDehYBhjDOFgmFAghDfdu/h8aV0R/pmFxa/JwWk8PjeldcWEAiHe+Md3AHB5nPiyfaRm+yjbU0zRtoLIa4YMTpetTjeqTWxTFQljzOSS20+LyDdEJN8YM2xlLrW5zIzOMto1ztTgNJODU0wNThMOGz7wx/cgDsE/s0A4bMivziUtNxVfVgq+LN/i/z/48J7rvn7h1vzrPp9Tnn3d53fcW7vykyIcfGQPcxPzzE3MMzsxx3S0aEGkq+u1f3ib1BwfmYXpZJVkklmcQVZxJi6P87r7VWo5m6pIiEgxMGCMMSJyhMjorBGLYykbW5gLMN47wXjPBFturcLtddF7cYCW19twepxkFqRTuruYzMJ0wmGD0yHsvO86H9IWc7ocFO8ovM7zTmpuq2J6ZIbx3kn6mgYBOPDwHkp2FjI7PsdIxxi5ldmkZvu0u0qtylZFQkS+B9wL5ItIN/BngBvAGPN3wKPAvxGRIDAHPGaM0a4kdY3Z8Tk6TnYz0j7K1NAMACJCQU0+OeVZVO4vpbSuKCE/JFMyvGy/u2bxvn9mgYn+SbJLswAYah2h8fnLAHjTveRWZpNXkU3JriJcXlt9HCibsNVPhTHmM6s8/9dEhsgqBUTOJ0wNzTDYMkxOeRZ5lTkEF4J0nuwhpzyLbXcVkVOedU13izfdi3eV100U3jTPNd1flYfKyK3KYbRrjNHOcUY6xuhvGqR4VxEAw+2jYCCnIlvPayjAZkVCqbUwxjDaNc5g8zADzUPMTcwDUHvnFvIqc8goSOeBf38XTrf2wb+biJCRn0ZGfhpVB8sxxjA/5ccdbUW0vtXOaOc4TreDvKpcirYXULQtH3eK29rgyjJaJNSmYMKG2fE50nJTATj7y0YWZgPkVeew9Wg1hbX5eNMiI4hERAvEGokIvsyUxfv1j+5npHOMoSsjDLYMM9gyTG91Dkc+HRneG1wI4vLox0Yy0e+2si1jDBP9U/Q1DtB3cQBjDPf94R04HA4OP7qP1GyffmBtMKfbSeHWfAq35lP3wHYm+qcw4chpP//sAi9/403yqnIo21tCYW0eTpcW40Snv2HKlgZbhrn4YjOzY3OIUyjcmk9JtN8cILMww8J0yUFEFicQAmCgqr6c3gsDDP38PC6vi9K6Impuq7qmNaISixYJZQvGGEbax/Blp5CWk4rT48Sb5qHm1iqKdxRon7gNeNM87Ly3lh13b2W4Y5Se8/30nO+j5tbIigaz43O4fe7F8xsqMeh3U1kqMB+g60wvnad6mJuYZ8uRSnbeV0teZQ55nztsdTy1DHEIBVvyKNiSRygQWjz/0/jcZUa7xynbXUzlwTIyCtItTqo2ghYJZZlLL7fQcbKHUCBEbkU22+/eStH2AqtjqXVYOkCg9s4tdJ7spvtsH52nesitzGbr0Wryq3MtTKhulhYJFVeTA1NkFkXOJ4QCYYp3FFBdX7H4mNq8sksyyf5IHTvvq6X7XB8dJ7uZ6JskvzqXcDgMBhxOnXux2WiRUDF3dV5DyxttjHaOc9u/OkxOWRa7HtiWcDOeFXhSI+eSqusrFkdG9TcNcunlVqpvqaRif6muI7WJaJFQMWOMYbh9lNY32xnrnoic+Ly/drGvWgtEYnM4HRCtBSmZKfiyfTS92MyVY+3U3FpF5cEync+yCWiRUDETCoQ589QFnG4ndQ9sp3xfiX4oJKnc8mxu++whxronaH79Ck0vtTDYMsytnz1kdTS1Ci0SakNNDU3TdbqXXe/bhsvj5JbHDpKel6brACkAcsqzOPLYQUY6xxa7ooILIQYuD1JaV4w4tHVpN1ok1IaYm5yn+bUr9Jzvx+V1UXGglIyCdLL0hLRaRl5lzuLtvosDnH+mibbjXey6fxt5VTnX+Z8q3rRIqJsSCoa5cqydK293goEtRyqpua0Kj08nv6m1udoNeemVFt75/ikKa/PZeV/t4jpdylpaJNRNEYG+i4MU1eaz475aXZ5BrZuIUFpXRNG2fNobumg91sGFZy8tXjNcWUuLhFq36eEZWt5sZ88Hd+Dyurj9d+v1gjXqpjndTrYeraZ8XylBfxAA/7SfsZ4JirYX6Gg4i+hvtlqzSNdSB61vteN0O5kamianPFsLhNpQ3jTP4rLv7Se6uXKsg4KaPOrev53UbN8q/1ttNP3tVmsy1jPB+V83MT0yQ8muIna9b9viL7JSsbLtri14Uj00v36F1/7H22y7cwvVt1TgcOhouXjRIqHWpPm1KwQXghx+dN81l8NUKpYcDgdbbqmgeEcBF5+/zKWXWwn6g2y/e6vV0ZKGFgm1orHuCXzZKaSke9n3kTpcHqd2LSlL+DJTOPixvfQ3DZIbHT7rn13Ak+LWuRUxpr/x6j3CoTAtb7TReqyD8r0l7P3QLlIyvFbHUklORBYvPGXChpM/OQvAvo/WkZajw2VjRTv21DWmh2d4659P0PpWB2V7Sth5/zarIyn1XgKVh8qZHpnljW8fp+dCv9WJEpa2JNSioSsjnPzpOZxuJwcf2UPxjkKrIym1LBGhbHcxuRXZnPnFBc7+spGR9lHq3r9dr3u+wfTdVIuySjIp2VnI9nu2kpKu3UvK/nyZKRz5zEFa32ynr2nQ6jgJSbubktzU8Axnf9VIOBTG43Oz7yN1WiDUpuJwONh2Zw13fPEWXB4XoUCIvosDVsdKGNqSSGI9F/q58JsmnG4ns2NzpOenWR1JqRvmdEWWoe881UPTSy0MR7ufrj6ubowWiSRkwoZLL7fQdryLnIpsDjy0W1sPKmFU11cQmA/Q+lYHU0PTHHxkr64pdhO0uykJXXjuMm3Hu6g6XM6RTx/QAqESijiE7Xdv5eDH9jI9Msub3znOWPe41bE2LW1JJKGqw+VkFWdQsb/U6ihKxUzx9gLS8+o5/+smPKm6hMyN0pZEkhjpGKXx+csYY8jIT9MCoZJCel4at37uEGm5qRhj6GsaxBhjdaxNRYtEEuht7Of4D88w0jFGcCFkdRyl4urqEuNDrSOc/vl5zjx1gVBAfw/WSotEAjPG0HqsgzO/aCSnPIvbPncIt669pJJUwdY8dty7lb6mQd75/ikW5gJWR9oUtEgksEsvt3L5lVZKdhVS/8kDuFP0kqIqeYkINbdWceDhPUwMTPH2v5xkfspvdSzb0z8rE1huZTYI7Lhnq17VS6mokp2FeHwuTv3sPNPDM7p45SokGU7i1NfXm4aGBqtjxEUoGGK0c5yCmjyroyhlawF/cLH7NTAfSLqWtoicMMbUr7addjclkOBCiBM/PsuJH59lZmzW6jhK2drVAjHQPMQrf/8WYz0TFieyJy0SCSLgD9Lwo9OMdI6x98M7dX19pdYosygDd4qbhh+e1kKxDC0SCSDgD3L8B6cZ753kwEN7KNtTYnUkpTaNqyvJelI9WiiWoUUiAfRfGmRyYIqDD++hZKdeA0Kp9VpaKI7/8DRzk/NWR7INHd2UACr2lZJTmqWruCp1E64Wiv6mQR3xtIStWhIi8i0RGRSR8ys8LyLy30SkRUTOisiheGe0i3AozJlfNjLRPwWgBUKpDeDLTGHLkUpEhKnhGaZHZqyOZDlbFQng28CD13n+Q8C26NfjwN/GIZPtmLDhzC8u0Huhn8nBKavjKJVwjDGc/tl5jv9Au55sVSSMMa8Co9fZ5GHgn0zEMSBbRJLqLK0xhosvNNN/aYid99dSsU8X6lNqo4kI+x+qI7gQ4vgPTuOfWbA6kmVsVSTWoAzoWnK/O/pY0rjydgcdJ7vZcqSSLbdUWh1HqYSVWZjB4Uf3MTc5T8OPzhD0B62OZInNViTWTEQeF5EGEWkYGhqyOs6GMMYw1j1BSV0RO+7danUcpRJebnk2hx7Zy9TgNK3HOqyOY4nNNrqpB6hYcr88+th7GGOeAJ6AyLIcsY8WeyLCoY/vBYOuxaRUnBRszaP+U/vJKc+yOoolNltL4ingd6OjnG4DJowxfVaHirXpkRne+cEp5qf8OBwOHM7N9m1TanPLr87F6XISmA/Qf2nQ6jhxZauWhIh8D7gXyBeRbuDPADeAMebvgKeBDwMtwCzwJWuSxk9gPsCJn5wl6A9iwmGr4yiV1FreaKe9oYtDH99L0bYCq+PEha2KhDHmM6s8b4A/jFMcy4XDYU797DxzE/Pc+pmD+LJ8VkdSKqltv7uGse5xzvyykaOfrycjCeYnab+FjTW92MJIxxh7PriDnPJsq+MolfScbieHPr4Pp8vBySfPEkiCEU9aJGwq4A8y3DZK9S0VlOtcCKVsIyXDy8FH9jA3MU/Ti81Wx4k5W3U3qd9ye13c/oV6HC6t40rZTW5FDgce2k1OWeKPeNJPIJsJ+oNcfrWVUCCEy+PC4dBvkVJ2VLyjEG+6FxM2zIwm7kW+9BPIRowxnH+midZjHUwOTFsdRym1Bheeu8yx755gftpvdZSY0CJhI52neuhrGmT7XTVJO3FHqc2m+nA5wYUQZ3/ZSGQAZmLRImETU0PTNL3YQkFNHjW3VVkdRym1Run5adQ9sJ2RjjGuJODSHVokbOLCs5dweZ3s/fAuXXJDqU2mfF8JJTsLaX6tLeEuf6qjm2xi30fqmJ+ax5vmsTqKUmqdRITdD+5kYS6AOBLrjzwtEhabn/LjTfeQmu0jNVtnVCu1Wbm9Lo48dtDqGBtOu5ssFPAHeeufG7j4QuJPyFEqWYQCIRqfv8xw2/Wun7Z5aJGw0OVXWpmf8lOyq8jqKEqpDTTcNsq5X19MiGU7tEhYZLRrjM5TPVTXVyTFrE2lkoXTHRmAMj/tT4hlO7RIWCAUCHHu1034slLYdleN1XGUUhsspyyLmiOVdJ/tY6RjzOo4N0WLhAVmx+cIBULseXAnLo/T6jhKqRiovWMLvuwULr7YvKkn2enoJgtkFKRzz+NHcbq1QCiVqJxuJ/s/uhuPz72p5z5pSyKOjDH0NvYTDoW1QCiVBHLKskjLTQUi3cybkRaJOOq7OMiZXzTS15Rc18hVKtmd/VUjJ588tym7nbRIxElwIUjTS81kFmVQqkNelUoqmUUZDLePMtg6YnWUddMiESetb7bjn16g7v3bE27avlLq+ioPlpGWm0rTi82EQ2Gr46yLFok4mB2fo+14F2V7inVOhFJJyOF0sPP+WmbH5ug42WN1nHXRIhEHwYUg2WVZbL97q9VRlFIWKajJI786l85T3Zjw5jk3oUNg4yCzMIPbPnvI6hhKKQuJCHs+FJkbtZm6nLUlEUPGGDpOdrMwu2B1FKWUDfgyU3CnuDHGENwk6zppkYih4bZRGp+7TO/FAaujKKVswoQNb/3zCRo3yerPWiRixBjD5Vdb8WWlULG/zOo4SimbEIeQU5ZFz/k+ZkZnrY6zKi0SMTLYOsLkwDS1d2zB6dK3WSn1WzW3VeF0OWh+vc3qKKvST68YMMbQ8kYbvuwUSut04pxS6lreNA9VhyvouzjA1NC01XGuS4tEDIQWQqRmpbD1aDUOp77FSqn32nKkEqfHSecpe8+b0CGwMeDyujj4yF6rYyilbMzjc3PbZw+Rnp9mdZTr0iKxwaaGphER23/jlVLWyyzKACJd1HZdTlz7QjZY00stvPODU4TDm2t9FqWUNYaujPDK37+Ff9pvdZRlaZHYQFND0wy3jVJ5sByHQ99apdTqUnN8zE3O097QZXWUZekn2QZqb+jC4XJQeaDU6ihKqU0iLSeV4h2FdJ7qIbhgv1nYWiQ2iH9mgd4LA5TtKcGT6rE6jlJqE6k+XEFwIURvo/1WZ9AisUEm+icRh1BdX251FKXUJpNdlklmUQYdJ7ptd/U6Hd20QQq35nP/V+/A5dG3VCm1PiLCjnvteSkB/UTbAMGFEC6PUwuEUuqG5VfnWh1hWdrdtAFO//w8DT8+Y3UMpdQmNz/t5+KLzcxP2Wc4rBaJmzQ3McfQlZHFSTFKKXWjQoEQ7ce76LnQb3WURVokblLX2T4AKvaVWJxEKbXZpeWkklOeRc+5PtucwLZVkRCRB0Xkkoi0iMjXlnn+iyIyJCKno19fsSLnVcYYes71kb8lF1+Wz8ooSqkEUba3hJnRWcZ7J62OAtioSIiIE/gb4ENAHfAZEalbZtMfGGMORL++GdeQ7zLaNc78lJ+yPcVWxlBKJZCSHYU43Q56zvVZHQWwUZEAjgAtxpgrxpgF4PvAwxZnuq6sogz2fWQXhbUFVkdRSiUIl9dF6e4SxGmPBf/sNGazDFi6eEk3cOsy231CRO4GLgN/bIxZdsETEXkceBygsrJyg6NGuLwuyvbouQil1Mba88EdVkdYZKeWxFr8Aqg2xuwDngO+s9KGxpgnjDH1xpj6goKN/0t/pHOMtnc6CQVCG/7aSikFkSGxVrNTkegBKpbcL48+tsgYM2KMufqufRM4HKds79F5socrb3fYpkmolEosza+38eoTb1n+h6idisRxYJuIbBERD/AY8NTSDURkad/OQ8DFOOZbFAqGGLoyQtH2Al0SXCkVE7kVWYQCYQZbRyzNYZtPOGNMEPgq8BsiH/4/NMZcEJE/F5GHopv9OxG5ICJngH8HfNGKrCMdY4QCIYq26QlrpVRs5FRk40l1039p0NIcdjpxjTHmaeDpdz32p0tu/yfgP8U717sNXB7C5XGSV5VjdRSlVIJyOBwUbS+g90I/oWAIp8tpTQ5L9rrJBRdCFNbm43Dq26eUip3C2nxCgTCjneOWZbBVS2KzOPjwHttMmVdKJa68yhz2/04d2aWZlmXQIrFOJmwQhyCio5qUUrHldDsprbN2RQftL1mnY989QePzl62OoZRKEguzC7S908ns+Jwl+9cisQ4LcwHGeyfx+NxWR1FKJYlgIETTSy0MNA9Zsn8tEusw0j4KQJ5NryCllEo8qVk+0nJTGekYs2T/WiTWYbh9FJfXRVaJXmBIKRU/uZXZjHWNEw6H475vLRLrMNI+Rl5Vjs6yVkrFVW5FDsGFEFMD03Hft45uWqNwOEzloTLS89OsjqKUSjK5ldmICNMjM2SVxHc4rBaJNXI4HNTcWmV1DKVUEkpJ9/LAf7gblyf+s66132SNpoam8c8sWB1DKZWkrCgQoEVizc4+fZHTT523OoZSKklN9E/yzvdPMTM6G9f9apFYg+BCkKmBaXLKsq2OopRKUg6ng5GOMcZ7J+K737jubZOa6J/CGEN2mXXrpyilklt6XhpOj5OJvqm47leLxBpMRoedZRVrkVBKWUMcQlZRBuN9k3HdrxaJNZganMKb5sGb5rE6ilIqiWWVZDI5OEUoGL9JdToEdg22Hq2mdLe1KzEqpVRuRTZTg9ME5gM4071x2acWiTVIy00lLTfV6hhKqSRXWJtPYW1+XPep3U2rmJuYo+tMLwtzAaujKKUUQFwveqZFYhWj3ROcf6aJBZ1Ip5SygRM/Ocupn8VvzpYWiVXMjM6CQGqOz+ooSimFw+VgcjB+w2C1SKxidnQWX1YKDqe+VUop66XnpzE3Pk8oEIrL/vSTbxUzY7Ok5ehJa6WUPaRFezXidTlTLRLXYYxhZnRORzYppWwjNTtaJCbm47I/HQJ7HSLCvX9wFBOO30gCpZS6ntScVMr3lcRtcq8WiVV4UnWWtVLKPjw+N3s/tCtu+9PupuuY6Juk+fUrOkdCKWUrxhgC8/H5XNIicR2j3RO0vNFudQyllLrGqZ+e4+3vnYrLvrRIXId/ah6Hy4E7RXvllFL24fa543alzHUXCRFJExFrrqMXZ/PTflIyvIiI1VGUUmqRN83LwuxCXAbVrFokRMQhIp8VkV+JyCDQBPSJSKOI/JWI1MY8pUXmp/x447TSolJKrZU3zQOGuJwvXUtL4iVgK/CfgGJjTIUxphC4EzgG/KWI/KsYZrRMYD6Ix+e2OoZSSl3DEx3+ujAb+y6ntXS2P2CMeU+5MsaMAj8BfiIiCflJeufvHSEcit/FPZRSai0yC9OpvXML7pTYf/SuWiSuFggRqQe+DlRF/59Enjb7lisiiUBEcLqS4vSLUmoTSctNZdsdW+Kyr/UM2/ku8L8C54CE//M6FAjR+PxlSuuKyavKsTqOUkotMsYwP+XH5XHGvDWxntFNQ8aYp4wxbcaYjqtfMUtmsYW5AN1n+5gZm7U6ilJKXSO4EOLlv32T7rN9Md/XeloSfyYi3wReAPxXHzTGPLnhqWwg6A8C4PbqHAmllL24PJFu8OBCMPb7Wse2XwJ2Am5+291kgIQsEuFg5BAdbj0noZSyFxHB6XYSXIj9NSXWUyRuMcbsiFkSm7k6qsmpFxtSStmQy+Nc7PGIpfV8Ar4pInUxS2Iz4bDB6XbicGuRUErZj8PliMuM6/W0JG4DTotIG5FzEotDYDcqjIg8CPx/gBP4pjHmL971vBf4J+AwMAJ82hjTvlH7XyqvMocP/Md7YvHSSil102rv2II3PfaXMlhPkXgwZimA6HpQfwO8H+gGjovIU8aYxiWbfRkYM8bUishjwF8Cn97oLMYYhq6MMDkwRWZRBgU1ebp+k1LKVsr3lsRlP+spEvnGmBNLHxCRjwIbNQz2CNBijLkSfe3vAw8DS4vEw8B/jt7+MfDXIiLGmA1rcxljOPnkOQZbhhcfK6zN59DH92qhUErZxszYLA6H4MvyxXQ/6+lw/wcR2XP1joh8Bvg/NjBLGdC15H539LFltzHGBIEJIG8DMzB0ZeSaAgEw2DLM0JWRjdyNUkrdlFM/Pc/FF5pjvp/1FIlHgX8SkZ0i8q+Bfwt8IDaxbp6IPC4iDSLSMDQ0tOb/NzkwtcLj0xsVTSmlbp7AxvWhrGzNRSLaDfQYkXkRnwA+YIyZ2MAsPUDFkvvl0ceW3UZEXEAWkRPYy+V9whhTb4ypLygoWHOIzKKMFR5PX/NrKKVUrJmwQRyx7wJf9ZyEiJwjMmnuqlwio4/eFhE2cHTTcWCbiGwhUgweAz77rm2eAr4AvEWkZfPiRp6PACioyaOwNv895yQKaja0V0sppW5KOBTGEYd5XGs5cf3RmKcgco5BRL4K/IZIEfqWMeaCiPw50GCMeQr4H8A/i0gLMEqkkGwoEeHQx/fS8mY7La+3seO+WrbcUqEnrZVStmJCBofTBi0JoHO1v9Y3aoSRMeZp4Ol3PfanS27PA5+82f2sRkTIKcsCILskQwuEUsp2dtxXG7lCXYyt6cp0IvJHIlK59EER8YjI/SLyHSJdQAklsyiDWz51gIwCPRehlLKfkp2F5FZkx3w/a2lJPAj8HvC96PmCcSCFSJfQs8B/Ncacil1Ea3h8bvK35FodQymlljXaNU5qto+UDG9M97OWK9PNA98AvhG9TGk+MGeMGY9pMouFgiEGW0bILEwnLTfV6jhKKbUoFAzz9r+cZNtdNdTeXh3Tfa3r1LgxJmCM6Uv0AgGRpcJP//w8g63Dq2+slFJxtHi9m5TYX+9GlzhdgcsTefPjsRSvUkqtR2A+AGiRsJQ4BKfHSdAf+4t6KKXUelz949UVhytnrrlIiMhzIrI/lmHsxp3iWqzYSillF4H5q5dXdsd8X+tpSfxvwH8VkX8UkfisUWsxb5qX+Wn/6hsqpVQcZRZlcPCRPaTlxX5QzZrbKsaYk8B9IvIJ4BkReRL4L8aYuZils9jeD+3Eqde4VkrZjDfNQ/GOwrjsa13nJCQy9fgS8LfAHwHNIvL5WASzg4yCdFKzY7tWu1JKrdd47yQjnWNx2dd6zkm8QWThvf+XyHUdvgjcCxwRkSdiEc5qU0MI9UxhAAAf1klEQVTTtB3vIhwKWx1FKaUWXXm7gwvPXorLvtZzavxxoHGZNZr+SEQubmAm2xjvmaDpxWaKdxTgy0yxOo5SSgHgn/aTkh7bmdZXred6Eheus4jfRzYoj62kRAvD3MS8xUmUUuq35qf9eO1WJK7n6nWpE01qTuR8xOx4wp6bV0ptMqFgmPkpP6nZ8end0Ml01+HLTEFEmB3TIqGUsoe58TkwkJoTnzXlYj9dbxNzOB2kZHqZHZ+1OopSSgHgy07h6Ofr8cWpJaFFYhW3fvYQntTYz2pUSqm1cLqcZJdmxm1/2t20Cl9mCk6XTqhTStnDQMswfRcH4rY/LRKrmBqapvH5y/hnFqyOopRSdDR00Xa8M2770yKxioW5AB0nupkcnLI6ilIqyRljmBycjutllbVIrOLqN2NqcMbiJEqpZOefWSAwF9AiYScenxtvuofp4WmroyilktzUYORzKKNQi4StZBSkMzWkRUIpZa3pkUiPRmYci4QOgV2DjIJ0BpqHMMYQWQhXKaXir7q+gpKdRbhT4jcsX4vEGmy/u4ad99VaHUMpleREhJSM+KzZdJV2N62Bw6lvk1LKWgtzAU4/dYGJvsm47lc//dao8fnLXHyx2eoYSqkkNd47Qd/FAYKBUFz3q0Vijean/Aw2D1sdQymVpMZ7JhERsorjtyQHaJFYs+zSTGbH5/DP6sxrpVT8jfdOkFGYjssT32WCtEisUU5ZNgBjXeMWJ1FKJZtwKMx432RcF/a7SovEGmWVZOB0O+N28XGllLrKP7NAapaPvKqcuO9bh8CukcPpoGxvMSkZeq1rpVR8+TJTuPP3jliyby0S67D7/TusjqCUSkJWTuTV7qZ1CofDBOYDVsdQSiWJ4EKIF/7763Sd6bVk/1ok1sEYw6tPHOPiCzpfQikVH2Pd4wTmAnGfaX2VFol1EBGyS7MYujKCCRur4yilksDQlREcTgc55dmW7F+LxDoV1uaxMBtgPM5T45VSyccYw2DLMHnVOXGfH3GVFol1KqjJQ0QYbNHZ10qp2JoenmFuYp7C2nzLMmiRWCd3ipuciiwtEkqpmHN5XWw9WkXhVuuKhA6BvQG1d2wBrB2WppRKfL7MFLbfvdXSDFokbkBeZfxnPSqlkot/ZoGJ/knyq3MtvVyBLbqbRCRXRJ4Tkebov8t+CotISEROR7+einfOpSYHp7jydoeVEZRSCay3cYATPz7L7PicpTlsUSSArwEvGGO2AS9E7y9nzhhzIPr1UPzivddw2xiXXm5lZmzWyhhKqQTVd3GAjMJ00vPSLM1hlyLxMPCd6O3vAI9YmGVNSnYVAtDfNGhxEqVUopkdn2Oib5LSuiKro9imSBQZY/qit/uBld6ZFBFpEJFjInLdQiIij0e3bRgaGtrQsBA5oZRdlkXfxYENf22lVHK7+rlSssv6IhG3E9ci8jxQvMxTX196xxhjRGSl6cxVxpgeEakBXhSRc8aY1uU2NMY8ATwBUF9fH5Pp0SW7Crn4fDOTg9NkFqbHYhdKqSQ02jVOTlkWvkzrV52OW5Ewxjyw0nMiMiAiJcaYPhEpAZbtwzHG9ET/vSIiLwMHgWWLRDyU7iqi9c12ZkZntUgopTZM/aP78c/Y4yqYduluegr4QvT2F4Cfv3sDEckREW/0dj5wB9AYt4TL8KR6uO8P76BkZ6GVMZRSCcQYgzjEsgX93s0uReIvgPeLSDPwQPQ+IlIvIt+MbrMLaBCRM8BLwF8YYywtEgAOhwNjjC4frpS6acGFEK/+wzFbneu0xWQ6Y8wI8L5lHm8AvhK9/SawN87R1qThR2fAwC2fPmB1FKXUJtZ/aZDZsTm8NmlFgH1aEptadmkmw+2jzE5YO+lFKbW5dZ/pJS03lZyyLKujLNIisQHK95UiInSe6rE6ilJqk5ocnGasZ2Lx88QutEhsAF9mCkXbC+g+00twIWR1HKXUJtTR0IXT7aBif4nVUa6hRWKDVNWXE5gP2uqEk1Jq8yjfX8qu923HneK2Oso1bHHiOhHklGVx6ON7KajJszqKUmoTyinLstW5iKu0JbFBRISibQWWLumrlNp8wqEwTS8223axUP1E22BdZ3o580vLp28opTaJ3sYB2o53MTtmz9GRWiQ2WGA+QO+Ffsb7Jq2OopSyORM2XDnWQUZhOvlbcq2OsywtEhus8kAZ7hQXrW+1Wx1FKWVz/ZcGmRmdZevRalsNe11Ki8QGc3ldVB0uZ7B5mKmhaavjKKVsyhhD67EO0nJTKd5eYHWcFWmRiIGqwxU43U5aj+nlTZVSywsHw2SXZlF7RzXisGcrAnQIbEx4fG523l9LWo7P6ihKKZtyup3s+eAOq2OsSotEjFQeKLM6glLKpkY6x3A4hJzybKujrEq7m2IoMB+g8fnLOtJJKbUoHA5z4dlLnHumCWNictHMDaVFIoZEhN7GAS6/YtnF85RSNtN7vp+ZkVm2373VtiOaltIiEUMur4utR6sY6RhjuH3U6jhKKYuFgiGaX28jqySTom35VsdZEy0SMVZ5sIyUDC+XXm7BhO3ftFRKxU7HiR7mp/zsuKdmU7QiQItEzDldTnbcu5XJgWm6z/VZHUcpZSGn20HJrkLyquw5u3o5OropDkp2FTE1NENupf1HMiilYqfqUDlVh8qtjrEu2pKIAxFhxz1bSctJtTqKUsoCE/1T9Fzo3xSjmd5Ni0QczU3Oc/LJs8yM2nNJYKXUxjPG0PjcJZpeaiG0Ca9cqUUijhwOYbhjjAvPXdqUf1Eopdav53w/472T7Li7Bpd38/Xwa5GII2+6lx33bGWkfYzeRr3MqVKJzj+7QNNLLWSXZVG2117Xrl4rLRJxVnmgjKySTC6+0MzCXMDqOEqpGGp6sYWgP8ieD+7YNENe302LRJyJQ9jz4E6C/iAtb7RZHUcpFUPF2wvYeV8tGQXpVke5YZuvgywBZBamc/CRPeRW6JBYpRJZkY2vE7FWWiQsUrQt8sMTCoYxYYPL47Q4kVJqo1x6uQWnJ7Isz2btZrpKu5ssFA6FOfY/G7j4wmWroyilNshIxyhX3u5kYWZh0xcI0CJhKYfTQX51Lt1n+xhsHbY6jlLqJgXmA5z91UXSclPZce9Wq+NsCC0SFqu9s4aMgjTO/bqJhdkFq+MopW7CxRea8U8vsO8jdTjdidGFrEXCYk6Xg30frSMwF+D8b3SSnVKb1fTIDL0XBth6tIrs0kyr42wYPXFtA5mFGWy/u4bexgGC/iDuFLfVkZRS65Sel8bRL9STkZ9mdZQNpUXCJrYcqaTqcDlOV2I0UZVKFuFQmLGeCfIqc8gqyrA6zobT7iabEBGcLicBf5CLLzYT9AetjqSUWoPLr7TyzvdOMTU0bXWUmNAiYTPTQ9O0N3Rx/lk9P6GU3Q20DNN2vIvKg2Wbelb19WiRsJmc8my23bmFvsYBOk/1WB1HKbWC2fE5zv2qkcyidHbeX2t1nJjRImFDW49WU7A1j4svNDPaNWZ1HKXUu4SCYU4+eRZj4MDDexL6XKIWCRsSEfb/zm5Ss31cePaydjspZTNOl4PKg+UceHh3wl9xUkc32ZTb6+LwJ/bhcDkSYmq/UoliYS6Ax+em8mCZ1VHiQlsSNpaWm4ovMwVjDH1Ng9qiUMpiA81DvPx3bzLWM2F1lLixRZEQkU+KyAURCYtI/XW2e1BELolIi4h8LZ4ZrdTfNMjpn5+n+XW9/oRSVpnon+TMLxpJz00lszAxRzItxxZFAjgPfBx4daUNRMQJ/A3wIaAO+IyI1MUnnrWKdxZSvq+E1jfb6Trba3UcpZLO7MQcDT8+i8fn5vAn9iXMukxrYYtzEsaYi8Bqfe9HgBZjzJXott8HHgYaYx7QYiLC7g/sYH7Sz4VnLpGS7qWgJs/qWEolhcB8gIYfniEcDHPrYwfxpnutjhRXdmlJrEUZ0LXkfnf0saTgcDo48Mge0gvSOPPLRoILOiNbqXhwepzkV+dy+BN7SU+wdZnWIm4tCRF5Hihe5qmvG2N+HoP9PQ48DlBZWbnRL28Jt9dF/aP7mR2fw+WxRSNQqYQVDocJzAXxpnmoe/92q+NYJm6fNMaYB27yJXqAiiX3y6OPrbS/J4AnAOrr6xNmWFBKhpeUjEhzt7dxgJzyLHyZKRanUiqxGGM493QTY93j3PGlI7i9yftH2WbqbjoObBORLSLiAR4DnrI4k2UW5gJcePYSx39wGv+MXqxIqY1ijKHxucv0XuinfF9pUhcIsEmREJGPiUg3cBT4lYj8Jvp4qYg8DWCMCQJfBX4DXAR+aIy5YFVmq3l8buof3cf81DzHf3CawHzA6khKbXrGGC6/0krnqR623FrJ1qNVVkeynCTDBK36+nrT0NBgdYyYGG4bpeEnZ8gsSKf+Uwfw+PSCRUrdqK4zvZx/ponKg2XUvX97Qq92ICInjDErzku7yhYtCXXj8rfkcuiRvUwOTTPYPGR1HKU2teIdBWy7qybhC8R6aJFIAIW1+dz9ldso31cKoMt3KLUOxhi6TvcQCoZwp7ipvb1aC8QSWiQSRGq2D4CJvkne/u5J/NN+ixMpZX/GGC6+0Mz531yi53y/1XFsSYtEggkGQkwOTnHsX04yOz5ndRylbCscDnP+mSY6TnRTXV9Bxf5SqyPZkhaJBJNXmcMtnzpIYC7Asf95gsnBKasjKWU7oUCIUz89T/fZPrYerWbn/bXaxbQCLRIJKKc8i1s/ewhEePu7JxP2Au1K3aj5aT/jvRPUvX872++u0QJxHck9SySBZRSkc/Tzh2l7p5O03MS+cpZSa7UwF8Cd4iItJ5W7Hz+a9BPl1kJbEgnMl5lC3QPbcTgd+Kf9dJzs1pFPKmlN9E/y+rfe5sqxDgAtEGukRSJJdJ7upfG5y5z7dRPhUNjqOErFVV/TIMe+exJxCAVb862Os6loKU0StXdUgzG0vNnO7Pgchx7ZgyfVY3UspWLKGEPrW+00v9ZGdmkmhz6+D2+a/tyvh7YkkoSIsO2uGvb/Th0TvZO8+c8NzIzOWh1LqZiaHp6h5fV2SncXceQzB7VA3ABtSSSZ0rpiUrN9NL3YgidV13lSiSkwH8Cd4o4M4PjdejKL0nUE0w3SlkQSyi7N4tbPHcKd4iYUDHPl7Q49T6ESxsDlIV7+u7cYbBkGIKs4QwvETdAikaSu/tIMtQ5z6eVWjn33JHOT8xanUurGhcNhml5u4eRPz5GW4yO9IPkuNRoLWiSSXPGOQg48vIfpkRne+PZxBnQlWbUJzU3O8873TtH2dicVB8q49XOHSc3yWR0rIWiRUJTsLOSOL9yCLzOFk0+eozU6jlypzWK4fZTJwWn2fbSOPR/cgdOlH20bRU9cKwDSclM5+vnDNL/eRmFtZBy5MUb7cpVtBfxBJgemyKvMoXxvCQU1eaSke62OlXC03KpFDqeDHfdsJSM/0pd77tdNtB7rIBzWk9rKXka7x3njH9/h5JPnCPiDiIgWiBjRloRaVigYJugPcvmVVgYuDbL3w7vIKEi3OpZKcsGFIJdfuULHyW58WSnUP7pPl9eIMX131bKcLgeHPraXvqZBGp+7xBvfPs7W26vZelsVDqc2QFX8BfxB3vjHd5ibmKfqcDnb767B5dGPsFjTd1hdV8nOQvIqs2l8vpmOhi4q95fi1Wa9iqNQMITT5cTtdVG+t4Tcqhxyy7OtjpU0tEioVXlSPRx4aDfzU3686V6MMVw51kHFgTI8Pp21rWLDGEP32T4uv9pK/Sf3k1WcSe0dW6yOlXS0SKg1S8mItCAm+iZpfq2NtuNd7LhnK+X7SnQUlNpQE/1TXHj2EhN9k+SUZ+FwOa2OlLS0SKh1yy7N4o4v3cKFZy9x/pkmus70sut928gpy7I6mkoATS+10PZOJ540D/s+sovS3cX6R4iF9AykuiEZBenc+tlD7PtoHfOT85x7+iImrBc0UjcmFAwtXhDL5XVSdbicu//1bZTt0Vaq1bQloW6YiFC2u5iibfnMTcwjDiG4EOTKsQ6qb6nU8xVqVSZs6LnQT/NrV6h7YDtF2wuovV3PO9iJFgl101we1+IcipHOcVqPddBxsoctt1RQVV+h49jVexhjGGwZpvm1NqaGpskqzsCj13qwJf3tVRuqqDafO790hMuvXqH59TbaG7qovqWSrbdVIQ7tNlARp352noHLQ/iyU9j/O7sp2VWo3Uo2pUVCbbiMgnQOf2IfE/2TNL/exmjnGLW3VwMQDoV1Ml4SMsYw2DpCflUOTreT0roiCmvzKa0r0p8Hm9MioWImqziT+kf3EwqGAJif8vP6P75D+Z5iquor8GWmWJxQxVooGKa3sZ/2411MD8+w58GdVOwvpXhHodXR1BppkVAx54yOcQ+Hw+RX59De0E37iW5KdhWx5UglmYW6JlSiCYfDtL3dSceJbvwzC2QUpLPvI7soqSuyOppaJy0SKm5Ss3wceGgPs/fM0X68i+6zffQ1DnD/V+/Ak6onLRPB/JSflAwvIsJgyzAZhensO1JJXlWOnnPYpLRIqLhLzfJR98B2tt25hZGOscUCceaXjfgyU6g4UKpdUZtIKBCir2mQrtM9TA5Mc98f3oHH5+bIYwdxunWm9GanRUJZxp3iXuybDociS5O3vtVO67F2CmryKNtTQmFt3mJ3lbKXucl52t7ppPdCP4H5IGm5qWy/pwZHdBSbFojEoEVC2YLD6eDwJ/YxOzFH1+lees73M9R6nl3v20Z1fQXhcBgR0S4Li81NzBEKhknPSyMcCtN5uofibQVUHCwjtyJbvz8JSK5OhU9k9fX1pqGhweoYah1M2DDSOUZGYTreVA/d5/poeaON4h2FFG0vILs0Uz+Q4mR+ys9A8xD9TYOMdo1TtL2AQx/bC0Su8aCTJTcnETlhjKlfbTv97ipbEoeQX527eD8lw0tabirtDV20vdOJN91L8fYCdr6vFodDx9nHyplfXKC3cQCIXAd9211bKK0rXnxeC0Ti0++w2hTyq3PJr84lMB9gsHWEgUuDTPRPLhaIK2934Pa5KdiSt7ikuVq7UCDEaNc4w22jjPVMcNvnDuFwOsgpzyYtL43i7QWkR699rpKLFgm1qbhT3JTtLqZsd/HiqqFXL04zMzoLQEZBGvlb8ijaXqDLl69itHucltfbGOueWJwNn1ORhX9mAV9mCpUHy6yOqCymRUJtWlfPSYgId33lVqaHZxi6MsJw2yjtJ7oiH3hlWQQXQlx6pYWcsmxyyrMWx/Enk+BCkPHeSca6xxnrnqDm1iryt0S68/zTC1QeLCN/Sy65Fdk6KkldQ4uESggiQkZBOhkF6dTcWkUoECIcCgMwMzZLz/l+Ok/2AOD2uckqzqD2ji3klGVFroMhJEzhCMwHCAXCpGR48U/7eft7pxZbWQhkFqYvLpWSW57NXV+51cK0yu5sUSRE5JPAfwZ2AUeMMcsORRKRdmAKCAHBtZyZV8nJ6XYu/kWcVZTBA//+LqYGpxnvnWSyf4qJganFotB/aZBzv24iLS+V9NxU0vLSSM9LJb86F5dNT8waYxbzd5zsZnp4hpnRWaZHZvFP+6nYX8qeB3fiSfWQUZBGya4isssyyS7N0pPNal3s8tNyHvg48Pdr2PY+Y8xwjPOoBONwOMgqziSrOPM9z/myUijfV8LMyCyjXeOLo3nu/Te34/K66DjRTdfZXnyZKfgyU/BmePD4PJTtKcbhdLAwFwAiV1S7mZFWxhhCCyEW5gKEQ5G5CAA9F/qZGpxmfmqeuUk/85PzZBSmU//ofgDa3ulcnMyWX51DWl4aeZXZQGSU2MFH9t5wJqVsUSSMMRchcZr7anPJLs0iu/S3J7iDC0FmRmcXR0l5Ut2kZHiZm5hjtGucoD8IQNneyFDQy69eoet0pCtLHILT7cTtc3Hv798OwKVXWhluG7lmn55UD7d86gAQubbCYOsw4WB48fn0vLTFbqCu0z2M903iy0ghJdNLXlXONcXuji8dweVx6u+PiglbFIl1MMCzImKAvzfGPLHShiLyOPA4QGVlZZziqUTg8riu+RAu2VVEya7frl4aXAgRmA8sthpK64pIz0sluBAkFAgTCoSueT13iouUjGvXolrajZVXnYMvKwWHy4HL48Tjc+NdMoy3/pP7cbpXLgLafaRiKW4zrkXkeaB4mae+boz5eXSbl4E/uc45iTJjTI+IFALPAX9kjHl1tX3rjGullLqW7WZcG2Me2IDX6In+OygiPwWOAKsWCaWUUjdm06xnICJpIpJx9TbwASInvJVSSsWILYqEiHxMRLqBo8CvROQ30cdLReTp6GZFwOsicgZ4B/iVMeYZaxIrpVRysMUZL2PMT4GfLvN4L/Dh6O0rwP44R1NKqaRmi5aEUkope9IioZRSakVaJJRSSq1Ii4RSSqkVaZFQSim1Ii0SSimlVqRFQiml1Iq0SCillFqRFgmllFIritsqsFYSkSGgYx3/JR/Y7Bc2SoRjgMQ4jkQ4BkiM40iEY4CNOY4qY0zBahslRZFYLxFp2OyXRk2EY4DEOI5EOAZIjONIhGOA+B6HdjcppZRakRYJpZRSK9IisbwVL4u6iSTCMUBiHEciHAMkxnEkwjFAHI9Dz0kopZRakbYklFJKrSjpi4SIfFJELohIWERWHC0gIu0ick5ETotIQzwzrsU6juNBEbkkIi0i8rV4ZlwLEckVkedEpDn6b84K24Wi34vTIvJUvHMuZ7X3VkS8IvKD6PNvi0h1/FNe3xqO4YsiMrTkvf+KFTlXIyLfEpFBEVn2EscS8d+ix3lWRA7FO+Nq1nAM94rIxJLvxZ/GJIgxJqm/gF3ADuBloP4627UD+VbnvZnjAJxAK1ADeIAzQJ3V2d+V8b8AX4ve/hrwlytsN2111vW+t8C/Bf4uevsx4AdW576BY/gi8NdWZ13DsdwNHALOr/D8h4FfAwLcBrxtdeYbOIZ7gV/GOkfStySMMReNMZesznGz1ngcR4AWY8wVY8wC8H3g4dinW5eHge9Eb38HeMTCLOuxlvd26bH9GHifiEgcM65mM/x8rIkx5lVg9DqbPAz8k4k4BmSLSEl80q3NGo4hLpK+SKyDAZ4VkRMi8rjVYW5QGdC15H539DE7KTLG9EVv9wNFK2yXIiINInJMROxQSNby3i5uY4wJAhNAXlzSrc1afz4+Ee2i+bGIVMQn2obbDL8La3FURM6IyK9FZHcsduCKxYvajYg8DxQv89TXjTE/X+PL3GmM6RGRQuA5EWmKVvq42aDjsNz1jmPpHWOMEZGVht9VRb8fNcCLInLOGNO60VnVe/wC+J4xxi8iv0+kZXS/xZmS1UkivwfTIvJh4GfAto3eSVIUCWPMAxvwGj3RfwdF5KdEmuZxLRIbcBw9wNK//Mqjj8XV9Y5DRAZEpMQY0xdt/g+u8BpXvx9XRORl4CCR/nSrrOW9vbpNt4i4gCxgJD7x1mTVYzDGLM37TSLnkDYjW/wu3AxjzOSS20+LyDdEJN8Ys6FrU2l30xqISJqIZFy9DXwAWHbEgc0dB7aJyBYR8RA5eWqLkUFLPAV8IXr7C8B7WkgikiMi3ujtfOAOoDFuCZe3lvd26bE9CrxoomcgbWLVY3hXv/1DwMU45ttITwG/Gx3ldBswsaSbc1MQkeKr57RE5AiRz/ON/6PD6jP4Vn8BHyPSH+kHBoDfRB8vBZ6O3q4hMtLjDHCBSPeO5dnXexzR+x8GLhP5q9uOx5EHvAA0A88DudHH64FvRm/fDpyLfj/OAV+2OvdK7y3w58BD0dspwI+AFuAdoMbqzDdwDP939HfgDPASsNPqzCscx/eAPiAQ/b34MvAHwB9Enxfgb6LHeY7rjGy08TF8dcn34hhweyxy6IxrpZRSK9LuJqWUUivSIqGUUmpFWiSUUkqtSIuEUkqpFWmRUEopG1ptgb93bXu3iJwUkaCIPPqu554RkXER+eWN5NAioZRS9vRt4ME1bttJZPHFf1nmub8CPn+jIbRIKKWUDZllFvgTka3RlsEJEXlNRHZGt203xpwFwsu8zgvA1I3m0CKhVAyIiE9EXhERp4hUr6XLYIXX8YjIq9FlPJR6AvgjY8xh4E+Ab8R6h/qDp1Rs/B7wpDEmdDOrgRtjFkTkBeDTwHc3KpzafEQknchqAz9a8jPljfV+tSWh1DqJyEsi8v7o7f9LRP77Mpt9juXXnaoRkVMicku0hdEkIt8Wkcsi8l0ReUBE3pDIlfmORP/bz6Kvp5KbAxg3xhxY8rUrHjtVSq3PnwFfF5HPEVl99j8sfTK6OF6NMab9XY/vAH4CfNEYczz6cC3w/wA7o1+fBe4k0pXwv0e3OQ/cEpMjUZuGiaz62iYin4TFS7Duj/V+tUgotU7RE4oC/EfgMWNM6F2b5APj73qsgEjL4nPGmDNLHm8zxpwzxoSJLNb2goksqHYOqI7uLwQsXF2JWCUHEfke8BawQ0S6ReTLRFqUXxaRq4uNPhzd9hYR6QY+Cfy9iFxY8jqvEVlY8n3R1/ngenLoOQml1klE9gIlwIgxZrlRI3NEVnxdaoLIMMU7uXZZc/+S2+El98Nc+/vpBeZvIrbaZIwxn1nhqfcMi422TMtXeJ27biaHtiSUWofo9RS+S+QvuGkRWe4XdgxwisjSQrFAZDn33xWRz65zn3nAsDEmcOPJlboxWiSUWiMRSQWeBP4XY8xF4P8kcn5iOc8SaTUsMsbMAB8F/lhEHlrHru8DfrX+xErdPL2ehFIxICKHgD82xtzwTNclr/Uk8DVjzOWbT6bU+mhLQqkYMMacBF4SEefNvE50pNTPtEAoq2hLQiml1Iq0JaGUUmpFWiSUUkqtSIuEUkqpFWmRUEoptSItEkoppVakRUIppdSKtEgopZRa0f8P1iY2mtSGJl8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from einsteinpy.plotting import StaticGeodesicPlotter\n",
    "sgp = StaticGeodesicPlotter(M)\n",
    "sgp.plot(pos_vec, vel_vec, end_lambda, stepsize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " All data regarding earth's orbit is taken from https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html"
   ]
  }
 ],
 "metadata": {
  "file_extension": ".py",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  },
  "mimetype": "text/x-python",
  "name": "python",
  "npconvert_exporter": "python",
  "pygments_lexer": "ipython3",
  "version": 3
 },
 "nbformat": 4,
 "nbformat_minor": 2
}